Epigenome-Wide Association Studies of the Fractional Exhaled Nitric Oxide and Bronchodilator Drug Response in Moderate-to-Severe Pediatric Asthma

Asthma is the most prevalent pediatric chronic disease. Bronchodilator drug response (BDR) and fractional exhaled nitric oxide (FeNO) are clinical biomarkers of asthma. Although DNA methylation (DNAm) contributes to asthma pathogenesis, the influence of DNAm on BDR and FeNO is scarcely investigated. This study aims to identify DNAm markers in whole blood associated either with BDR or FeNO in pediatric asthma. We analyzed 121 samples from children with moderate-to-severe asthma. The association of genome-wide DNAm with BDR and FeNO has been assessed using regression models, adjusting for age, sex, ancestry, and tissue heterogeneity. Cross-tissue validation was assessed in 50 nasal samples. Differentially methylated regions (DMRs) and enrichment in traits and biological pathways were assessed. A false discovery rate (FDR) < 0.1 and a genome-wide significance threshold of p < 9 × 10−8 were used to control for false-positive results. The CpG cg12835256 (PLA2G12A) was genome-wide associated with FeNO in blood samples (coefficient= −0.015, p = 2.53 × 10−9) and nominally associated in nasal samples (coefficient = −0.015, p = 0.045). Additionally, three CpGs were suggestively associated with BDR (FDR < 0.1). We identified 12 and four DMRs associated with FeNO and BDR (FDR < 0.05), respectively. An enrichment in allergic and inflammatory processes, smoking, and aging was observed. We reported novel associations of DNAm markers associated with BDR and FeNO enriched in asthma-related processes.


Introduction
Asthma is a heterogeneous respiratory disease that affects more than 330 million people worldwide and is the most prevalent chronic disease in children. Moreover, asthma is associated with one in every 250 deaths worldwide [1]. Characterized by chronic airway inflammation and variable expiratory airflow limitation, asthma is defined by its symptoms, including wheezing, shortness of breath, chest tightness, and cough, which may vary over time and in intensity. Anti-inflammatory and bronchodilator drugs are the main therapy to control asthma symptoms, improve lung function, and reduce the risk of severe exacerbations [2]. However, less than 5% of pediatric asthma patients have severe asthma, defined as uncontrolled asthma that often does not respond to the currently available medications despite adherence to maximal optimized high doses [3]. Severe asthma patients have the lowest quality of life and the highest risk for morbidity and mortality, and their treatment consumes the majority of healthcare resources due to asthma [3].
The fractional exhaled nitric oxide (FeNO) and the bronchodilator drug response (BDR) are two of the main clinical biomarkers that can be objectively measured for asthma diagnosis and for assessing treatment response. FeNO is the only available non-invasive biomarker for type 2 airway inflammation in asthma. It can be quantified in exhaled air of children and adults through portable and non-invasive devices. FeNO levels in exhaled air have been related to chronic inflammatory diseases, Th2-mediated immune responses, and airway eosinophilia [2]. BDR measures changes in lung function after administration of short-acting β 2 -agonists [4], the main rapid-onset airway bronchodilators. BDR is also associated with the severity and instability of asthma and is also frequently used for asthma phenotyping. FeNO and BDR are conditioned by genetic and environmental factors [5,6].
Precision medicine aims to improve the prevention and treatment of diseases based on features inherent to each patient [7]. The genetic background greatly contributes to asthma, estimating the heritability of asthma susceptibility over 55-90% and 70% to treatment response [7]. However, currently reported genetic variants only partially explain this heritability. Asthma is a multifactorial disease, and susceptibility and treatment response are influenced by genetic variation and environmental factors [7]. DNA methylation (DNAm), an epigenetic mechanism that reflects the interaction between both genetic and environmental factors, could explain part of the heritability of asthma treatment response [8]. Epigenome-wide association studies (EWAS) have identified promising CpG sites and differentially methylated regions (DMRs) involved in asthma [9,10]. Despite being useful clinical biomarkers in the management of asthma, the relationship of DNAm with FeNO and BDR has been scarcely investigated. Cardenas et al. [9] profiled DNAm in nasal cells in 547 children from different populations and identified a large set of CpGs and DMRs associated with multiple asthma-related phenotypes, including FeNO (8,372 CpGs and 191 DMRs with FDR < 0.05) and BDR (130 CpGs with an FDR < 0.05). They reported that FeNOassociated epigenetic markers were enriched in intracellular membrane trafficking, T cell activation, oxidative stress, mucin production, neutrophil degranulation, and interleukin pathways. Multiple DMRs were located in asthma-associated genes, including TNIP-1, IL-13, and CHI3L1. For the BDR trait, they did not report any significant DMRs and observed that the associated epigenetic markers did not overlap with any other asthmarelated phenotype. However, this study only included 12% of individuals with current asthma. On the other hand, FeNO-related genes, including multiple nitric oxide synthase (NOS) and arginase (ARG) isoforms, have been evaluated in candidate-gene epigenetic studies, but controversial results of the effect of epigenetic markers in these genes and FeNO have been reported [11,12].
The association of FeNO and BDR with both environmental and genetic factors has been reported, but its relationship with the epigenome in asthma patients remains barely investigated. Thus, we hypothesize that DNAm from whole blood is a biomarker associated with BDR and FeNO that captures the influence of environmental and genetic factors on these clinical biomarkers. Therefore, this study aimed to identify DNAm markers in whole blood associated with BDR or FeNO in pediatric asthma, assess whether they are shared with nasal epithelium cells, and analyze their enrichment in previous associations and/or biological processes.

Study Population
This study was conducted in children with moderate-to-severe asthma from the Systems Pharmacology Approach to Uncontrolled Pediatric Asthma (SysPharmPediA) consortium (Clinicaltrials.gov ID: NCT04865575) [13]. This study was conducted following the Declaration of Helsinki and approved by the ethics committees of participating institutions: University Regensburg, Germany (18- Medical Ethics Committee, Slovenia (0120-569/2017/4); and Ethical Regional Review Board in Stockholm, Sweden (Dnr 2010/1336-31/3 and 2019-00546). All parents and participants provided signed written consent and agreed to participate in the study. The study design and characteristics of the enrolled patients have already been described elsewhere [13]. Briefly, patients were enrolled in four European countries (Spain, The Netherlands, Slovenia, and Germany) according to the following inclusion criteria: (1) aged between 6 and 17 years old, (2) physician's diagnosis of asthma, and (3) moderate-to-severe patients with asthma treated with medication step ≥ 3 according to the Global Initiative for Asthma (GINA) guidelines [2]. In the current study, only individuals with available data on BDR and/or FeNO were included.
Clinical and demographic variables of interest were recorded using standardized questionnaires. Asthma control was assessed based on (childhood) Asthma Control Test (ACT/cACT) and the development of asthma exacerbations. Uncontrolled asthma was defined as having an ACT/cACT score ≤ 19 or ≥ 1 exacerbation requiring oral corticosteroids or severe exacerbation requiring hospitalization or emergency room visits in the past year. Controlled asthma was defined as asthma patients with an ACT/cACT score > 19 and no severe exacerbations in the past year. Body mass index (BMI) was estimated as BMI = weight (kg)/height (m 2 ). According to the World Health Organization (WHO) recommendations, we expressed BMI as z-scores [14].

Measurement of FeNO and BDR
Spirometry tests pre-and post-administration of bronchodilator and FeNO assessment were performed according to the European Respiratory Society/American Thoracic Society (ERS/ATS) guidelines [15,16], as previously described [13]. BDR was calculated based on the forced expiratory volume in the first second (FEV 1 ), expressed in liters, as BDR = postFEV1−preFEV1 preFEV1 × 100. FeNO values were reported as parts per billion (ppb).

Genome-Wide Methylation Assessment and Quality Control
Blood samples were collected in K-EDTA vacutainer tubes for subsequent downstream applications, including genotyping and methylation profiling. Genomic DNA was isolated using the FlexiGene DNA Kit (Qiagen) according to the manufacturer's instructions. After concentration measurements and quality control (QC), the samples were equilibrated to 10 ng/µL, and a total of 500 ng DNA was utilized for bisulfite conversion. Next, methylation levels of 865,918 CpG sites across the whole genome were assessed using the Infinium Illumina MethylationEPIC BeadChip array (Illumina Inc., San Diego, CA, USA). Raw data QC was conducted using the ENmix package (1.28.8) [17] in R (4.1.2) [18]. First, we removed low-quality samples and probes. Methylation data points were defined as low-quality data if their signal did not differ from negative control (detection p-value > 1 × 10 −6 ) and/or they were measured by < 3 beads. Based on that, low-quality CpG probes were defined as those having ≥ 5% of bad-quality data points across samples. Low-quality samples were defined as (1) the presence of ≥ 5% of bad quality data points across CpGs, (2) having total bisulfite intensity less than three standard deviations than bisulfite controls, and/or (3) outliers of bisulfite intensity or beta value distribution.
We used the out-of-band (oob) method for background noise correction. Bisulfite intensities were normalized by quantile normalization to reduce variability among arrays. The EPIC array uses two different dyes and types of probes; (i) Regression on Correlated Probes (RCP) and (ii) the Regression on Logarithm of Internal Control (RELIC) methods [19,20] to correct dye-and probe-type biases. Methylation intensities were used to estimate beta values. Furthermore, outlier methylation data points for each CpG, defined as values below/above three times the interquartile range (IQR) from the first/third quantiles, were set as missing values. CpGs and samples with a missing rate of ≥5% and ≥10%, respectively, were removed. The remaining missing values were imputed (k-nearest neighbor method). Moreover, we estimated the predicted sex based on methylation data and discarded one sample with sex discordance between the reported and predicted sex.
We removed first-and second-degree related individuals previously identified based on genotype data obtained with the Global Screening Array (GSA, Illumina Inc.). Since we included both males and females, we removed probes within sex chromosomes. Potential cross-sample contamination was assessed based on the distribution of beta-values of the control genotype probes included in the EPIC array. We used the ewastools R package [21] to identify potential cross-sample contaminated individuals and inspected the beta value distribution plots. Potential problematic probes that may capture other artifacts different than methylation were also filtered out. These include probes: (1) with a multimodal distribution of beta values, (2) cross-reactive or non-specific, which bind non-specifically to the target region of interest [22], and (3) potentially polymorphic probes defined as those containing a single nucleotide polymorphism (SNP) with a minor allele frequency (MAF) > 1% based on the Illumina v1.0 B4 manifest file at the CpG site or a single base extension. Finally, beta values were transformed into M-values for improved statistical properties [23].

Estimation of Cell-Type Heterogeneity and Potential Confounders
Cell-type heterogeneity is the major confounder in EWAS performed in primary tissues. In whole blood, DNAm captures the epigenetic states of a mixture of cell types. These distinct methylation profiles may lead to false discoveries when looking for their correlation with the phenotype of interest [24]. Therefore, cell-type heterogeneity in whole blood was captured using the Reference-Free Adjustment for Cell-Type composition (ReFACTor) [24] algorithm in GLINT 1.0.4, adjusting for sex, age, and the first six genotype principal components (PCs). This method estimates a set of PCs correlated with cell-type composition and could be included as covariates in subsequent analyses to correct this major confounder. A principal component regression analysis was used to test for association between multiple potential confounding variables and global DNAm.

Epigenome-Wide Association Study
We conducted an EWAS to identify single CpG sites across the genome associated with FeNO or BDR. The association between the methylation M-values and FeNO/BDR was tested through linear regression models using the limma v3.48.3 R package [25], adjusting for age, sex, ancestry (six genotype PCs), and the first six ReFACTor components. A false discovery rate (FDR) of 10% was used to correct for multiple comparisons. The p-value < 9 × 10 −8 threshold was used to declare genome-wide significant association, as recommended for the EPIC array [26]. Genome inflation factor (λ) and quantile-quantile plots (Q-Q plots) of regression p-values were examined to inspect genomic inflation. Probes were annotated to the nearest gene based on the Illumina manifest file v1.0 B4 [27] and GREAT v4.0.4 [28]. Furthermore, we evaluated the robustness of the results by conducting sensitivity analyses adjusting for ICS dosage and recruiting center and stratified analysis by sex, age, and ethnicity. Potentially problematic probes were filtered out for subsequent analyses.

Cross-Tissue Validation
We attempted to evaluate the effect of the CpGs identified in whole blood samples in 62 nasal samples. Clinical data are summarized in Supplementary Table S1. Nasal epithelial cells were collected using sterile nasal swabs (Paul Böttger OHG, Bodenmais, Germany) and stored at −80 • C until processing. Genomic DNA was isolated using the MasterPure Complete DNA & RNA Purification Kit (Lucigen, Teddington, UK) according to the manufacturer's instructions. Genome-wide DNAm profiling was carried out using the Illumina Infinium MethylationEPIC array. QC of DNAm data, cell-tissue heterogeneity assessment, and association analyses were conducted following the same procedures described for the blood samples.

Differentially Methylated Regions
A DMR analysis was conducted to identify genomic regions of several CpG sites in which the global DNAm might be associated with our outcomes of interest. Two independent software, comb-p [29] and DMRcate [30], were used to identify DMRs associated with BDR and FeNO. Regions were identified based on CpGs with p-value < 0.05 (comb-p) or FDR ≤ 0.2 (DMRcate), and those CpGs located within 750 bp (comb-p) or 1000 bp (DMRcate) were combined into the same region. Only DMRs significant after correcting for multiple testing overlapping in both methods (adjusted p-value < 0.05) were retained. Comb-p was used to identify DMRs associated with BDR and FeNO in nasal samples using the same parameters as in the discovery phase.

Enrichment Analyses
Enrichment analysis was next performed to identify whether a set of CpGs or genes were more likely to be related to biological pathways, gene ontology terms, or different traits than by chance. The top 100 probes from the EWAS of FeNO and BDR were included in CpG-set enrichment analyses to identify an enrichment of CpGs previously associated with any traits through EWAS. These analyses were carried out using the EWAS Toolkit [31]. Furthermore, the genes to which these CpGs were annotated were included in gene-set enrichment analyses to identify enrichment in biological pathways and gene ontologies using Enrichr [32].

Quality Control of DNAm Data and Assessment of Global DNAm Patterns
After QC, epigenomic data were available for 124 individuals and 773,260 CpGs from whole blood, and 55 individuals and 636,555 CpGs from nasal epithelia were available for subsequent analyses. The individuals and CpGs removed or flagged during the QC are summarized in Supplementary Table S2. The principal regression analysis plot showed that age, sex, ReFACTor principal components (capturing cell-heterogeneity), and ancestry had the most significant associations with global methylation patterns ( Supplementary Figures S1 and S2). As a result, these variables were considered the most relevant confounders to be included as covariates in the association analyses. FeNO and BDR also showed significant associations with global methylation patterns.

Epigenome-Wide Association Study
We did not observe genomic inflation both in the EWAS of BDR (λ = 1.04) ( Figure 1A) and FeNO (λ = 1.03) ( Figure 1C), indicating an optimal control for cell heterogeneity and other confounders. A total of three CpGs were significantly associated with BDR with FDR < 0.1 (Table 2, Figure 1B). The top hit was the CpG cg26203256, located in the ADD3-AS1 gene, in which DNAm was associated with lower BDR (coefficient = −0.02, p = 1.85 × 10 −7 ). On the other hand, the CpG cg12835256 located near PLA2G12A was genome-wide significantly associated with a reduction in FeNO (coefficient = −0.015, p = 2.53 × 10 −9 ) (Table 2, Figure 1D). The negative coefficients for both associations suggest that these CpG are more likely to be hypomethylated in patients with higher FeNO or BDR values. These associations remained robust in all sensitivity analyses considering ICS dosage, recruiting centers, and ethnicity (Supplementary Table S3). The top-hit CpG cg12835256 (PLA2G12A) was cross-tissue associated in nasal samples (coefficient = −0.015, p = 0.045, Supplementary Table S4). This CpG is flagged in Illumina's manifest file as a potential polymorphic probe. However, the potentially implicated SNP (rs4557260) is almost monomorphic in our population (minor allele count = 1). associations remained robust in all sensitivity analyses considering ICS dosage, recruiting centers, and ethnicity (Supplementary Table S3). The top-hit CpG cg12835256 (PLA2G12A) was cross-tissue associated in nasal samples (coefficient = −0.015, p = 0.045, Supplementary  Table S4). This CpG is flagged in Illumina's manifest file as a potential polymorphic probe. However, the potentially implicated SNP (rs4557260) is almost monomorphic in our population (minor allele count = 1).  Given that the age of the patients recruited in the study included puberty onset (13 years), where asthma prevalence trend switches between males and females [33,34], we also performed sensitivity analyses stratifying by sex and age. In the sex-stratified analysis, with the exception of the CpG from C21orf91, the rest of the CpGs were significant and showed consistent effects in both sexes (Supplementary Table S5). In the age-stratified anal- ysis, all CpGs showed consistent effects in both age groups, but the one from FAM69B was not statistically significant in patients younger than 13 years old (Supplementary Table S5).

Discussion
This study investigated the role of DNAm in two clinical biomarkers, FeNO and BDR, in children with moderate-to-severe asthma. We identified one CpG located in the PLA2G12A gene as genome-wide significantly associated with FeNO and three CpGs significantly associated with BDR. Furthermore, twelve and four DMRs were related to FeNO and BDR, respectively. In addition, the top-hit CpG and DMR associated with FeNO and BDR, respectively, were cross-tissue validated in nasal samples. CpG-set and gene-set analyses revealed significant enrichment in previous EWAS signals and biological pathways involved in asthma. Our findings are relevant by providing new perspectives on the influence of DNAm on asthma biomarkers in pediatric asthma. Although some results overlap with a previous study (e.g., biological pathways such as interleukin signaling), our CpGs and differentially methylated sites and also its implications in biological pathways had not been previously reported by Cardenas et al. [9].
The CpG genome-wide, significantly associated with FeNO in whole blood and crossvalidated in nasal samples (cg12835256), is located at the gene encoding for phospholipase A2 Group XIIA (PLA2G12A), a member of a group of phospholipases with a key role as regulators of type-2 inflammation, airway hyperresponsiveness, and eicosanoids production in asthma [35]. This enzyme, highly expressed in airway cells [36], acts as a high-affinity ligand for the PLA2R1 receptor, whose expression is elevated in children with asthma compared to healthy patients [37]. We also reported another CpG associated with FeNO, annotated to the chromosome 21 open reading frame 91 (C21orf91) gene. In a previous study, this gene was observed to be co-expressed with the IL-6 receptor (IL-R6), which has been associated with asthma risk. [38]. The complex IL-6/IL6R triggers a cascade that activates STAT3 and CEBP transcription factors, both involved in airway inflammation and asthma [39]. Increased levels of IL-6 have been related to impaired altered lung function and the use of high-dose ICS [40].
The most significant DMR for FeNO was annotated to PPP5D1, a pseudogene whose role in asthma is still unclear. Furthermore, we reported a DMR associated with FeNO located near EGR3. This gene encodes for a protein member of the early growth factors (EGR), which regulate the expression of genes related to inflammation and cell growth [41]. In addition, EGR3 was reported to be overexpressed in the lung of asthma patients and regulates the expression of IL-2, IL-6, and IL-8 [42][43][44]. The other DMRs associated with FeNO are located in genes of the HOXA gene family (i.e., HOXA-AS3, HOXA7, and HOXA4), AURKC, MAGI2-AS3, PDZRN4, and MOGAT3. HOXA-AS3 encodes a long noncoding RNA identified as a positive regulator of the necrosis factor kappa B (NF-κB) signaling [45]. NF-κB is a well-known transcription factor that leads to the expression of multiple pro-inflammatory mediators (e.g., cytokines and chemokines) involved in asthma [46]. In addition, MOGAT3 encodes a protein that catalyzes the synthesis of diacylglycerol, a metabolite that attenuates the pathological consequences of inflammatory processes [47]. Furthermore, DNAm markers in HOXA7 and HOXA4 genes have been associated with lung development and asthma [48], while DNAm in AURKC has been associated with lung function in Latino children with asthma [49]. Finally, MAGI2-AS3 encodes a lncRNA previously associated with asthma-related functions and pathways, such as cytokine production, cell proliferation, and kinase activity [50].
We reported the association of three epigenetic markers with BDR, located in ADD3-AS1, PARP2, and FAM69B. The ADD3-AS1 gene is located in the antisense strand of ADD3, whose gene expression patterns correlate with serum immunoglobulin E (IgE) levels in atopic asthma patients [51]. The PARP2 encodes a protein from PARP enzymes, which play key roles in asthma pathogenesis by affecting the expression of pro-inflammatory genes and chemokines [52]. Furthermore, we reported DMRs associated with BDR annotated to C5orf63 (validated in nasal samples), PPFIBP2, LY6G5C, and RAMP1. The receptor activity-modifying protein 1 (RAMP1) is a component of the Calcitonin gene-related protein (CGRP) receptor. Studies questioning the role of CGRP during allergic asthma have been contradictory, suggesting both protective (e.g., maintenance of vascular tone) and inflammation-promoting roles [53].
The epigenetic markers identified in this study showed enrichment in traits and biological processes, including PPARs, interleukins (IL-5 and IL-2), FcERI signaling pathways, smoking, and aging. PPARs have been established as critical components of type 2 immune response to allergens. The mechanism of action of PPARs is characterized by the ability to regulate lipid metabolism and inhibit proinflammatory transcription factors (e.g., NF-κB), which could make them a therapeutic target for asthma [54]. Patients with severe asthma have abnormal production of type 2 cytokines, whereas IL-5 is the main mediator of eosinophilic inflammation [55]. Similarly, alterations in the FcERI pathway are described in allergic diseases and associated with asthma exacerbations and allergic inflammation [56]. Interestingly, the strongest enrichment signal was related to aging. In the context of asthma, epigenetic age acceleration has been associated with asthma, allergic phenotypes, IgE, and FeNO [9,57]. In addition, we identified enrichment in CpGs related to smoking and second-hand smoke, which are well-known risk factors for asthma, impaired lung function, and unresponsiveness to corticosteroids [58].
We acknowledge some strengths of this work. Our study on pediatric patients provides a novel insight into a specific asthma phenotype that has been widely understudied: moderate-to-severe asthma. Moreover, we have controlled our analyses for known confounders in EWAS, such as age, sex, tissue heterogeneity, and ancestry. We also used microarrays that allow the agnostic genome-wide screening of DNAm changes with a low technical variation [59]. In addition, the multi-center pan-European nature makes the results likely generalizable to multiple populations. However, some limitations of this study must be acknowledged. First, the sample size of our study is modest in comparison to previous EWAS, not allowing us to detect DNAm changes with small effects. However, the prevalence of children with moderate to severe asthma is generally low, making this a unique cohort. Second, focusing on this specific phenotype of moderate-to-severe pediatric asthma could limit the generalizability of the results to other asthma phenotypes. Third, we have not included an independent population to replicate our findings, but we assessed cross-tissue validation in nasal samples. Fourth, although we assessed DNAm in whole blood and nasal samples, other asthma-relevant tissues remained unexplored. Fifth, although we adopted the standard approach of adjusting by cell heterogeneity, we acknowledge that this could have reduced the statistical power to detect DNAm changes related to eosinophilic or Th2 asthma. Sixth, although we addressed any potential bias in our results related to age and sex, we have not recorded sex hormone levels or the age of puberty onset to better classify the patients in sensitive and stratified analyses.

Conclusions
In conclusion, we reported novel associations of epigenetic markers in whole blood associated with FeNO and BDR in children with moderate-to-severe asthma. These markers were enriched in asthma-related pathways, such as aging, smoking, and inflammatory and allergic responses. Furthermore, these findings provide new insights into potential epigenetic biomarkers of asthma that could be used as predictors of prognosis or treatment response.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biomedicines11030676/s1, Table S1: Clinical and demographic characteristics of the individuals with methylation measured in nasal samples; Table S2: Summary table of the quality control of methylation data; Table S3: Sensitivity analysis for the main EWAS findings; Table S4: Evaluation of the effect of CpGs in nasal epithelial cells; Table S5: Summary results of sensitivity analyses stratifying by sex and age; Table S6: Summary results of the gene-set enrichment analysis for the EWAS of BDR; Figure S1: Principal component regression analysis plot of blood samples; Figure S2: Principal component regression analysis plot of nasal samples; Figure  S3: CpG-set enrichment analysis plot for the top 100 CpGs associated with BDR. Supplementary SI: SysPharmPediA consortium members (authors list); Supplementary SII: Acknowledgements.